#source('http://bioconductor.org/biocLite.R')
#biocLite('phyloseq')
library(phyloseq)
library(ggplot2)
library(dplyr)
library(Rmisc)
setwd("~/Dropbox/clado-manuscript/Mikes_MS_Data/")
# Load biom file. 
biom <- import_biom("OTU_table.biom", "~/Dropbox/clado-manuscript/Nephele/PipelineResults_NMEPINZ20QK1/nephele_outputs/tree.tre", parseFunction=parse_taxonomy_greengenes)
sam.data <- read.csv(file="sample.data.csv", row.names=1, header=TRUE)
head(sam.data)
##        TreatmentGroup  Site Date                       Description
## C172N1          Early North  172 Sample of day 172 at site North 1
## C172N2          Early North  172 Sample of day 172 at site North 2
## C172N3          Early North  172 Sample of day 172 at site North 3
## C172P1          Early Point  172 Sample of day 172 at site Point 1
## C172P2          Early Point  172 Sample of day 172 at site Point 2
## C172P3          Early Point  172 Sample of day 172 at site Point 3
##        SampleID.1
## C172N1     C172N1
## C172N2     C172N2
## C172N3     C172N3
## C172P1     C172P1
## C172P2     C172P2
## C172P3     C172P3
sam.data$Date <- as.factor(sam.data$Date)
sam.data$DateSite <- paste(sam.data$Date, sam.data$Site)
head(sam.data); str(sam.data)
##        TreatmentGroup  Site Date                       Description
## C172N1          Early North  172 Sample of day 172 at site North 1
## C172N2          Early North  172 Sample of day 172 at site North 2
## C172N3          Early North  172 Sample of day 172 at site North 3
## C172P1          Early Point  172 Sample of day 172 at site Point 1
## C172P2          Early Point  172 Sample of day 172 at site Point 2
## C172P3          Early Point  172 Sample of day 172 at site Point 3
##        SampleID.1  DateSite
## C172N1     C172N1 172 North
## C172N2     C172N2 172 North
## C172N3     C172N3 172 North
## C172P1     C172P1 172 Point
## C172P2     C172P2 172 Point
## C172P3     C172P3 172 Point
## 'data.frame':    52 obs. of  6 variables:
##  $ TreatmentGroup: Factor w/ 2 levels "Early","Late": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Site          : Factor w/ 3 levels "North","Point",..: 1 1 1 2 2 2 3 3 3 1 ...
##  $ Date          : Factor w/ 6 levels "172","178","185",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ Description   : Factor w/ 52 levels "Sample of day 172 at site North 1",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ SampleID.1    : Factor w/ 52 levels "C172N1","C172N2",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ DateSite      : chr  "172 North" "172 North" "172 North" "172 Point" ...
sample_data(biom) <- sam.data
biom; sample_data(biom)
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 51928 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 51928 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 51928 tips and 51914 internal nodes ]
## Sample Data:        [52 samples by 6 sample variables]:
##        TreatmentGroup  Site Date                       Description
## C178N1          Early North  178 Sample of day 178 at site North 1
## C178P1          Early Point  178 Sample of day 178 at site Point 1
## C185P2          Early Point  185 Sample of day 185 at site Point 2
## C206N2           Late North  206 Sample of day 206 at site North 2
## C206P1           Late Point  206 Sample of day 206 at site Point 1
## C206P2           Late Point  206 Sample of day 206 at site Point 2
## C214P1           Late Point  214 Sample of day 214 at site Point 1
## C214P2           Late Point  214 Sample of day 214 at site Point 2
## C214P3           Late Point  214 Sample of day 214 at site Point 3
## C214S1           Late South  214 Sample of day 214 at site South 1
## C214S2           Late South  214 Sample of day 214 at site South 2
## C214S3           Late South  214 Sample of day 214 at site South 3
## C185P1          Early Point  185 Sample of day 185 at site Point 1
## C185P3          Early Point  185 Sample of day 185 at site Point 3
## C199P3           Late Point  199 Sample of day 199 at site Point 3
## C199S2           Late South  199 Sample of day 199 at site South 2
## C199S3           Late South  199 Sample of day 199 at site South 3
## C206N1           Late North  206 Sample of day 206 at site North 1
## C178P2          Early Point  178 Sample of day 178 at site Point 2
## C199N3           Late North  199 Sample of day 199 at site North 3
## C206S1           Late South  206 Sample of day 206 at site South 1
## C214N3           Late North  214 Sample of day 214 at site North 3
## C199N2           Late North  199 Sample of day 199 at site North 2
## C206N3           Late North  206 Sample of day 206 at site North 3
## C206S2           Late South  206 Sample of day 206 at site South 2
## C199N1           Late North  199 Sample of day 199 at site North 1
## C199P1           Late Point  199 Sample of day 199 at site Point 1
## C199S1           Late South  199 Sample of day 199 at site South 1
## C214N1           Late North  214 Sample of day 214 at site North 1
## C172P1          Early Point  172 Sample of day 172 at site Point 1
## C199P2           Late Point  199 Sample of day 199 at site Point 2
## C172N3          Early North  172 Sample of day 172 at site North 3
## C172S3          Early South  172 Sample of day 172 at site South 3
## C178S2          Early South  178 Sample of day 178 at site South 2
## C178P3          Early Point  178 Sample of day 178 at site Point 3
## C178S3          Early South  178 Sample of day 178 at site South 3
## C172N1          Early North  172 Sample of day 172 at site North 1
## C172S1          Early South  172 Sample of day 172 at site South 1
## C178N3          Early North  178 Sample of day 178 at site North 3
## C185N2          Early North  185 Sample of day 185 at site North 2
## C185N3          Early North  185 Sample of day 185 at site North 3
## C185S3          Early South  185 Sample of day 185 at site South 3
## C214N2           Late North  214 Sample of day 214 at site North 2
## C172P2          Early Point  172 Sample of day 172 at site Point 2
## C185S2          Early South  185 Sample of day 185 at site South 2
## C172P3          Early Point  172 Sample of day 172 at site Point 3
## C185N1          Early North  185 Sample of day 185 at site North 1
## C172N2          Early North  172 Sample of day 172 at site North 2
## C178S1          Early South  178 Sample of day 178 at site South 1
## C185S1          Early South  185 Sample of day 185 at site South 1
## C172S2          Early South  172 Sample of day 172 at site South 2
## C178N2          Early North  178 Sample of day 178 at site North 2
##        SampleID.1  DateSite
## C178N1     C178N1 178 North
## C178P1     C178P1 178 Point
## C185P2     C185P2 185 Point
## C206N2     C206N2 206 North
## C206P1     C206P1 206 Point
## C206P2     C206P2 206 Point
## C214P1     C214P1 214 Point
## C214P2     C214P2 214 Point
## C214P3     C214P3 214 Point
## C214S1     C214S1 214 South
## C214S2     C214S2 214 South
## C214S3     C214S3 214 South
## C185P1     C185P1 185 Point
## C185P3     C185P3 185 Point
## C199P3     C199P3 199 Point
## C199S2     C199S2 199 South
## C199S3     C199S3 199 South
## C206N1     C206N1 206 North
## C178P2     C178P2 178 Point
## C199N3     C199N3 199 North
## C206S1     C206S1 206 South
## C214N3     C214N3 214 North
## C199N2     C199N2 199 North
## C206N3     C206N3 206 North
## C206S2     C206S2 206 South
## C199N1     C199N1 199 North
## C199P1     C199P1 199 Point
## C199S1     C199S1 199 South
## C214N1     C214N1 214 North
## C172P1     C172P1 172 Point
## C199P2     C199P2 199 Point
## C172N3     C172N3 172 North
## C172S3     C172S3 172 South
## C178S2     C178S2 178 South
## C178P3     C178P3 178 Point
## C178S3     C178S3 178 South
## C172N1     C172N1 172 North
## C172S1     C172S1 172 South
## C178N3     C178N3 178 North
## C185N2     C185N2 185 North
## C185N3     C185N3 185 North
## C185S3     C185S3 185 South
## C214N2     C214N2 214 North
## C172P2     C172P2 172 Point
## C185S2     C185S2 185 South
## C172P3     C172P3 172 Point
## C185N1     C185N1 185 North
## C172N2     C172N2 172 North
## C178S1     C178S1 178 South
## C185S1     C185S1 185 South
## C172S2     C172S2 172 South
## C178N2     C178N2 178 North
head(otu_table(biom))
## OTU Table:          [6 taxa and 52 samples]
##                      taxa are rows
##                                C178N1 C178P1 C185P2 C206N2 C206P1 C206P2
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      1
## New.CleanUp.ReferenceOTU321320      0      0      0      0      0      0
## KC551585.1.1451                     2      0      0      6      0      4
## JQ945994.1.1399                     0      0      0      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                   137      8      1     73      8      6
##                                C214P1 C214P2 C214P3 C214S1 C214S2 C214S3
## New.CleanUp.ReferenceOTU155901      0      0      1      0      1      0
## New.CleanUp.ReferenceOTU321320      0      0      0      0      0      0
## KC551585.1.1451                     9      4     11      2      1      0
## JQ945994.1.1399                     1      0      0      1      0      0
## EF653577.1.1339                     0      0      0      4      0      0
## JQ814729.1.1495                    11      4     27     40      8     23
##                                C185P1 C185P3 C199P3 C199S2 C199S3 C206N1
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      1      0      0      0      0
## KC551585.1.1451                     3      2      5      0      0     11
## JQ945994.1.1399                     0      0      0      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                     6      4      2     11      2     56
##                                C178P2 C199N3 C206S1 C214N3 C199N2 C206N3
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      0      0      0      0      0
## KC551585.1.1451                     0      1      4      3      0      3
## JQ945994.1.1399                     0      0      2      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                    10     23      6     19      5    849
##                                C206S2 C199N1 C199P1 C199S1 C214N1 C172P1
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      0      0      0      0      0
## KC551585.1.1451                     1      0      8      2      5      0
## JQ945994.1.1399                     0      0      0      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                    14     24      5      5      3      2
##                                C199P2 C172N3 C172S3 C178S2 C178P3 C178S3
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      0      0      0      0      0
## KC551585.1.1451                     2      5      8      0      1      0
## JQ945994.1.1399                     0      0      0      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                     2    121     31      0      0      5
##                                C172N1 C172S1 C178N3 C185N2 C185N3 C185S3
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      0      0      0      0      0
## KC551585.1.1451                     0      0      1      0      0      0
## JQ945994.1.1399                     0      0      0      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                     5     16      4      6     27      4
##                                C214N2 C172P2 C185S2 C172P3 C185N1 C172N2
## New.CleanUp.ReferenceOTU155901      0      0      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      0      1      0      0      0
## KC551585.1.1451                     3      2      0      2      1      6
## JQ945994.1.1399                     0      0      0      0      0      0
## EF653577.1.1339                     0      0      0      0      0      0
## JQ814729.1.1495                     4     13      1      3     22     21
##                                C178S1 C185S1 C172S2 C178N2
## New.CleanUp.ReferenceOTU155901      0      0      0      0
## New.CleanUp.ReferenceOTU321320      0      0      0      0
## KC551585.1.1451                     0      1      0      3
## JQ945994.1.1399                     0      0      3      0
## EF653577.1.1339                     0      0      0      0
## JQ814729.1.1495                     2      3     21     18
# Custom plotting. 
nolegend <- theme(legend.position="none")
readabund <- labs(y="read abundance")
# Normalize by relative abundance. 
biom.relabund <- transform_sample_counts(biom, function(x) x / sum(x))
ordNMDS <- ordinate(biom.relabund, method="NMDS", distance="bray")
## Run 0 stress 0.1164511 
## Run 1 stress 0.1164499 
## ... New best solution
## ... Procrustes: rmse 0.0002640801  max resid 0.001355786 
## ... Similar to previous best
## Run 2 stress 0.1164499 
## ... Procrustes: rmse 4.59707e-06  max resid 2.549379e-05 
## ... Similar to previous best
## Run 3 stress 0.1164499 
## ... Procrustes: rmse 2.024418e-06  max resid 7.244113e-06 
## ... Similar to previous best
## Run 4 stress 0.1164499 
## ... New best solution
## ... Procrustes: rmse 4.512229e-07  max resid 1.37136e-06 
## ... Similar to previous best
## Run 5 stress 0.1164511 
## ... Procrustes: rmse 0.0002640564  max resid 0.001356755 
## ... Similar to previous best
## Run 6 stress 0.1164511 
## ... Procrustes: rmse 0.0002640336  max resid 0.001355855 
## ... Similar to previous best
## Run 7 stress 0.1164511 
## ... Procrustes: rmse 0.000264165  max resid 0.001352902 
## ... Similar to previous best
## Run 8 stress 0.2490019 
## Run 9 stress 0.1164499 
## ... Procrustes: rmse 1.262105e-05  max resid 3.529697e-05 
## ... Similar to previous best
## Run 10 stress 0.1164499 
## ... Procrustes: rmse 3.669909e-06  max resid 1.715039e-05 
## ... Similar to previous best
## Run 11 stress 0.1164499 
## ... Procrustes: rmse 2.022824e-06  max resid 1.086017e-05 
## ... Similar to previous best
## Run 12 stress 0.1164499 
## ... Procrustes: rmse 8.11061e-07  max resid 2.358948e-06 
## ... Similar to previous best
## Run 13 stress 0.1164499 
## ... Procrustes: rmse 6.679075e-07  max resid 1.942404e-06 
## ... Similar to previous best
## Run 14 stress 0.1164499 
## ... Procrustes: rmse 2.77383e-06  max resid 1.271664e-05 
## ... Similar to previous best
## Run 15 stress 0.1164511 
## ... Procrustes: rmse 0.0002640178  max resid 0.001354858 
## ... Similar to previous best
## Run 16 stress 0.1164511 
## ... Procrustes: rmse 0.0002638534  max resid 0.00135605 
## ... Similar to previous best
## Run 17 stress 0.1164499 
## ... Procrustes: rmse 6.057614e-06  max resid 3.703111e-05 
## ... Similar to previous best
## Run 18 stress 0.1164511 
## ... Procrustes: rmse 0.0002642032  max resid 0.00135728 
## ... Similar to previous best
## Run 19 stress 0.1164511 
## ... Procrustes: rmse 0.0002639396  max resid 0.001359735 
## ... Similar to previous best
## Run 20 stress 0.1164499 
## ... Procrustes: rmse 2.275713e-06  max resid 9.427331e-06 
## ... Similar to previous best
## *** Solution reached
ordNMDS.k3 <- ordinate(biom.relabund, method="NMDS", distance="bray", k=3)
## Run 0 stress 0.07869986 
## Run 1 stress 0.07870007 
## ... Procrustes: rmse 0.0002212905  max resid 0.0006546112 
## ... Similar to previous best
## Run 2 stress 0.07869804 
## ... New best solution
## ... Procrustes: rmse 0.0004127187  max resid 0.002047968 
## ... Similar to previous best
## Run 3 stress 0.08345599 
## Run 4 stress 0.07869827 
## ... Procrustes: rmse 0.0002309752  max resid 0.001112312 
## ... Similar to previous best
## Run 5 stress 0.08696054 
## Run 6 stress 0.07869976 
## ... Procrustes: rmse 0.0004640888  max resid 0.001623611 
## ... Similar to previous best
## Run 7 stress 0.07890412 
## ... Procrustes: rmse 0.00682703  max resid 0.03194959 
## Run 8 stress 0.07869838 
## ... Procrustes: rmse 0.0003182201  max resid 0.001253299 
## ... Similar to previous best
## Run 9 stress 0.08337992 
## Run 10 stress 0.08337912 
## Run 11 stress 0.07899677 
## ... Procrustes: rmse 0.008309483  max resid 0.03642162 
## Run 12 stress 0.07870018 
## ... Procrustes: rmse 0.0006233716  max resid 0.002629041 
## ... Similar to previous best
## Run 13 stress 0.07875932 
## ... Procrustes: rmse 0.003019153  max resid 0.01347439 
## Run 14 stress 0.08337967 
## Run 15 stress 0.07886131 
## ... Procrustes: rmse 0.005660624  max resid 0.02758483 
## Run 16 stress 0.08695834 
## Run 17 stress 0.07869963 
## ... Procrustes: rmse 0.0005062642  max resid 0.001462566 
## ... Similar to previous best
## Run 18 stress 0.07870052 
## ... Procrustes: rmse 0.0006595409  max resid 0.002172851 
## ... Similar to previous best
## Run 19 stress 0.07869934 
## ... Procrustes: rmse 0.0003273543  max resid 0.001333105 
## ... Similar to previous best
## Run 20 stress 0.08695869 
## *** Solution reached
ord.k3 <- plot_ordination(biom.relabund, ordNMDS.k3, shape="Site", color = "DateSite") + geom_point(size=2) + geom_polygon(aes(fill=DateSite), alpha=0.7) + labs(title = "Cladophora, 2014") + theme_bw()
ord.k3

# Facet by Date. 
ord.k3.facet1 <- plot_ordination(biom.relabund, ordNMDS.k3, shape="Site", label = "Date") + geom_point(size=2.5) + facet_wrap(~Date) + labs(title = "Cladophora, 2014") + geom_polygon(aes(fill=DateSite)) + theme_bw()
ord.k3.facet1

# Facet by Site. 
ord.k3.facet2 <- plot_ordination(biom.relabund, ordNMDS.k3, label = "Date") + geom_point(size=2.5) + facet_wrap(~Site, ncol = 1) + labs(title = "Cladophora, 2014") + geom_polygon(aes(fill=DateSite)) + theme_bw()
ord.k3.facet2

# ANOSIM...
# Remove singleton. (EDA)
biom.nosingle <- prune_taxa(taxa_sums(biom)>1, biom)
biom.nosingle # Same, so QIIME QC covered it. 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 51928 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 51928 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 51928 tips and 51914 internal nodes ]
# Find methanotrophs
methanolist <- read.table(file = "~/Dropbox/clado-manuscript/clado_16S-archive/methanos.txt")
methanolist <- as.vector(methanolist$V1)
# 
biom.relabund.methanos <- subset_taxa(biom.relabund, Genus %in% as.factor(methanolist))
biom.relabund.methanos
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 567 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 567 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 567 tips and 566 internal nodes ]
head(tax_table(biom.relabund.methanos))
## Taxonomy Table:     [6 taxa by 8 taxonomic ranks]:
##                                Kingdom    Phylum          
## AB240506.1.1496                "Bacteria" "Proteobacteria"
## KC432108.1.1352                "Bacteria" "Proteobacteria"
## New.CleanUp.ReferenceOTU247624 "Bacteria" "Proteobacteria"
## AJ422152.1.1406                "Bacteria" "Proteobacteria"
## AB989868.1.1453                "Bacteria" "Proteobacteria"
## AJ422161.1.1420                "Bacteria" "Proteobacteria"
##                                Class                Order            
## AB240506.1.1496                "Betaproteobacteria" "Burkholderiales"
## KC432108.1.1352                "Betaproteobacteria" "Burkholderiales"
## New.CleanUp.ReferenceOTU247624 "Betaproteobacteria" "Burkholderiales"
## AJ422152.1.1406                "Betaproteobacteria" "Burkholderiales"
## AB989868.1.1453                "Betaproteobacteria" "Burkholderiales"
## AJ422161.1.1420                "Betaproteobacteria" "Burkholderiales"
##                                Family           Genus         Species
## AB240506.1.1496                "Comamonadaceae" "Methylibium" NA     
## KC432108.1.1352                "Comamonadaceae" "Methylibium" NA     
## New.CleanUp.ReferenceOTU247624 "Comamonadaceae" "Methylibium" NA     
## AJ422152.1.1406                "Comamonadaceae" "Methylibium" NA     
## AB989868.1.1453                "Comamonadaceae" "Methylibium" NA     
## AJ422161.1.1420                "Comamonadaceae" "Methylibium" NA     
##                                Rank1
## AB240506.1.1496                NA   
## KC432108.1.1352                NA   
## New.CleanUp.ReferenceOTU247624 NA   
## AJ422152.1.1406                NA   
## AB989868.1.1453                NA   
## AJ422161.1.1420                NA
# 
relabund.methanos <- psmelt(biom.relabund.methanos)
relabund.methanos.genus <- relabund.methanos%>%
  group_by(Sample, Genus)%>%
  mutate(GenusAbundance = sum(Abundance))%>%
  distinct(Sample, GenusAbundance, TreatmentGroup, Site, Date, Family, Genus)
head(relabund.methanos.genus)
## Source: local data frame [6 x 9]
## Groups: Sample, Genus [6]
## 
##   Sample GenusAbundance TreatmentGroup   Site   Date          Family
##    <chr>          <dbl>         <fctr> <fctr> <fctr>          <fctr>
## 1 C178S2       1.369438          Early  South    178 Crenotrichaceae
## 2 C199S1       1.369438           Late  South    199 Crenotrichaceae
## 3 C185S1       1.369438          Early  South    185 Crenotrichaceae
## 4 C199S3       1.369438           Late  South    199 Crenotrichaceae
## 5 C185S2       1.369438          Early  South    185 Crenotrichaceae
## 6 C178S1       1.369438          Early  South    178 Crenotrichaceae
## # ... with 3 more variables: Genus <fctr>, Sample <chr>, Genus <fctr>
# 
p <- ggplot(relabund.methanos.genus, aes(as.factor(Date), GenusAbundance, color = Site))
p <- p + geom_point() + facet_wrap(~Genus, scales="free_y")
p

# Means with error bars
# stats <- summarySE(biom, measurevar="Abundance", groupvars=c("bacsp","media")); stats
# p.stats <- ggplot(stats, aes(x = media, y = percloss, fill = bacsp))
# p.stats + geom_bar(stat = "identity", position=position_dodge(.9)) + geom_errorbar(aes(ymin=percloss-se, ymax=percloss+se), width=.2, colour="darkblue", position=position_dodge(.9)) + geom_rug()  + scale_fill_grey()

# Plot richness. 
biom.rich <- plot_richness(biom, x="Date", color="Site")
biom.rich

### Shannon linear model test. 
#biom.rich[1][1]

# Stacked bar plots of methanotrophs. 
barstack.methanos <- plot_bar(biom.relabund.methanos, x = "Date", fill="Site") + facet_wrap(~Genus, scales="free_y")
barstack.methanos

# Network plot. 
network <- plot_net(biom, maxdist = 0.3, point_label = "SampleID.1", color = "Date", shape = "Site")
network

# Trees :)
biom.relabund.Methylococcaceae <- subset_taxa(biom.relabund, Family %in% "Methylococcaceae")
sample_data(biom.relabund.Methylococcaceae)$Date <- as.factor(sample_data(biom.relabund.Methylococcaceae)$Date)
biom.relabund.Methylococcaceae
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 79 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 79 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 79 tips and 78 internal nodes ]
head(tax_table(biom.relabund.Methylococcaceae))
## Taxonomy Table:     [6 taxa by 8 taxonomic ranks]:
##                                Kingdom    Phylum          
## New.CleanUp.ReferenceOTU16149  "Bacteria" "Proteobacteria"
## New.CleanUp.ReferenceOTU412315 "Bacteria" "Proteobacteria"
## New.CleanUp.ReferenceOTU12803  "Bacteria" "Proteobacteria"
## New.CleanUp.ReferenceOTU418181 "Bacteria" "Proteobacteria"
## New.CleanUp.ReferenceOTU42851  "Bacteria" "Proteobacteria"
## New.CleanUp.ReferenceOTU410894 "Bacteria" "Proteobacteria"
##                                Class                 Order            
## New.CleanUp.ReferenceOTU16149  "Gammaproteobacteria" "Methylococcales"
## New.CleanUp.ReferenceOTU412315 "Gammaproteobacteria" "Methylococcales"
## New.CleanUp.ReferenceOTU12803  "Gammaproteobacteria" "Methylococcales"
## New.CleanUp.ReferenceOTU418181 "Gammaproteobacteria" "Methylococcales"
## New.CleanUp.ReferenceOTU42851  "Gammaproteobacteria" "Methylococcales"
## New.CleanUp.ReferenceOTU410894 "Gammaproteobacteria" "Methylococcales"
##                                Family             Genus             
## New.CleanUp.ReferenceOTU16149  "Methylococcaceae" "Methylomonas"    
## New.CleanUp.ReferenceOTU412315 "Methylococcaceae" "Methylomicrobium"
## New.CleanUp.ReferenceOTU12803  "Methylococcaceae" "Methylomicrobium"
## New.CleanUp.ReferenceOTU418181 "Methylococcaceae" "Methylomonas"    
## New.CleanUp.ReferenceOTU42851  "Methylococcaceae" "Methylomonas"    
## New.CleanUp.ReferenceOTU410894 "Methylococcaceae" "Methylomonas"    
##                                Species      Rank1
## New.CleanUp.ReferenceOTU16149  NA           NA   
## New.CleanUp.ReferenceOTU412315 "buryatense" NA   
## New.CleanUp.ReferenceOTU12803  NA           NA   
## New.CleanUp.ReferenceOTU418181 NA           NA   
## New.CleanUp.ReferenceOTU42851  NA           NA   
## New.CleanUp.ReferenceOTU410894 NA           NA
tree1 <- plot_tree(biom.relabund.Methylococcaceae, color="Date", shape="Date", label.tips="Genus", size = "Abundance")
tree1

tree2 <- plot_tree(biom.relabund.Methylococcaceae, color="DateSite", label.tips = "Genus", size = "Abundance")
tree2

# Sandbox.